load "localsplit.txt";

//given a number, n, and a prime, p, 
//returns the result of dividing all factors of p from n 
//followed by the highest power of p diving n
function extract_p(n,p)
	n := Integers()!n;
	p := Integers()!p;
	return (n div (p^Valuation(n,p))),Valuation(n,p);
end function;


function parse_matrix(M)
	N := GetLocal(M,2);
	Ei := 1;
	Mi := 1;
	Ni := 1;
	e_set := [];
	m_set := [];
	n_set := [];
	el_set := [];
	ml_set := [];
	nl_set := [];
	skip := 0;
//Break up the GetLocal matrix output into the three parts
//diagonalized part are the e's
//[0 1/2, 1/2 0] are the M's
//[1 1/2, 1/2 1] are the N's
	for i in [1..NumberOfRows(N)] do
		if (skip eq 1) then
			skip := 0;
		else
			if (i eq NumberOfColumns(N)) then
				e_set[Ei],el_set[Ei] := extract_p(N[i,i],2);
				Ei := Ei + 1;
			elif (N[i,i] eq 0) then
				m_set[Mi], ml_set[Mi] := extract_p(2*N[i+1,i],2);
				Mi := Mi + 1;
				skip := 1;
			elif (N[i+1,i] ne 0) then
				n_set[Ni], nl_set[Ni] := extract_p(2*N[i+1,i],2);
				Ni := Ni + 1;
				skip := 1;
			else
				e_set[Ei],el_set[Ei] := extract_p(N[i,i],2);
				Ei := Ei + 1;
			end if;
		end if;
	end for;
	return e_set, el_set, m_set, ml_set, n_set, nl_set;
end function;


function L(el_set,k)
	L_set := [];
	Li := 1;
	e_indices := [];
	for i in [1..#el_set] do
		if (el_set[i]-k lt 0) then
			if (IsOdd(el_set[i]-k)) then
				L_set[Li] := el_set[i];
				e_indices[Li] := i;
				Li := Li + 1;
			end if;
		end if;
	end for;
	return Li-1, L_set, e_indices;
end function;


function p(nl_set,k)
	p_ans := 1;
	for i in [1..#nl_set] do
		if (nl_set[i] lt k) then
			p_ans := p_ans * (-1)^(nl_set[i] - k);
		end if;
	end for;
	return p_ans;
end function;


function e(e_set,el_set,k)
	Li, L_set, e_indices := L(el_set,k-1);
	e_ans := 1;
	for i in [1..#e_indices] do
		e_ans := e_ans * e_set[e_indices[i]];
	end for;
	return e_ans;
end function;


function d(el_set, ml_set, nl_set, k)
	el_sum := 0;
	ml_sum := 0;
	nl_sum := 0;
	for i in [1..#el_set] do
		if (el_set[i] lt k - 1) then
			el_sum := el_sum + el_set[i]-k+1;
		end if;
	end for;
	for i in [1..#ml_set] do
		if (ml_set[i] lt k) then
			ml_sum := ml_sum + ml_set[i]-k;
		end if;
	end for;
	for i in [1..#nl_set] do
		if (nl_set[i] lt k) then
			nl_sum := nl_sum + nl_set[i]-k;
		end if;
	end for;
	d_ans := k + 1/2 * el_sum + ml_sum + nl_sum;
	return d_ans;
end function;


function delta(el_set,k)
	for i in [1..#el_set] do
		if (el_set[i] eq k-1) then
			return 0;
		end if;
	end for;
	return 1;
end function;

function hilbert(a)
	a := Integers()!a;
	if (a mod 8) eq 1 then return 1;
	end if;
	if (a mod 8) eq 7 then return 1;
	end if;
	if (a mod 8) eq 3 then return -1;
	end if;
	if (a mod 8) eq 5 then return -1;
	end if;
	return 0;
end function;

function mu(k,alpha,a,el_set,e_set)
	sum := 0;
	for i in [1..#el_set] do
		if el_set[i] lt k-1 then
			sum := sum + e_set[i];
		end if;
	end for;
	return alpha*2^(a-k+3)-sum;
end function;

function char_(mu)
	if (mu mod 4) eq 0 then
		return 1;
	end if;
	return 0;
end function;

function psi(n)
// We only care when char_(mu) is non-zero, so we assume mu = 0 mod 4
//	return Exp(-2*Pi(ComplexField())*Sqrt(-1)*n);

	if (Integers()!(8*n) mod 8 eq 4) then 
		return -1;
	end if;
	return 1;
end function;

function local_density_2(M,t)
	e_set, el_set, m_set, ml_set, n_set, nl_set := parse_matrix(M);
	alpha, a := extract_p(t,2);
	density := 0;
	for k in [1..a+3] do
		mu_ := mu(k,alpha,a,el_set,e_set);
		if (L(el_set,k-1) mod 2) eq 1 then
			density := density + delta(el_set,k)*p(nl_set,k)*hilbert(mu_*e(e_set,el_set,k))*2^(Integers()!(d(el_set, ml_set, nl_set, k)-3/2));
		else
			density := density + delta(el_set,k)*p(nl_set,k)*hilbert(e(e_set,el_set,k))*2^(Integers()!(d(el_set, ml_set, nl_set, k)-1))*psi(mu_/8)*char_(mu_);
		end if;
	end for;
	return 1+density;
end function;